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ABSTRACT 

Given the locations of several landmarks on a satellite acquired image and their true geo- 
graphic coordinates, the position and orientation of the satellite can be determined. Two 
methods for automatically locating the image coordinates of specified landmarks are des- 
cribed. The first, a particularly fast sequential similarity detection algorithm for template 
matching was originally described by Nagel and Rosenfeld. The second method involves 
iteratively resampling the picture function in the vicinity of the anticipated landmark. A 
variety of other speedup methods is also described. An application to SMS imagery is 
envisioned. 

INTRODUCTION 

An important role in satellite navigation is played by the recognition and location of land- 
marks on images of the earth’s surface as viewed by the satellite. Given the locations of 
several landmarks on the image, and their true geographical locations, the position and 
orientation of the satellite can be determined. The mathematics of this process, and its 
error analysis, will not be recapitulated here. (See, for example, Phillips and Smith, 1972.) 

We shall assume here that the approximate position and orientation of the satellite are 
already known, and that we want to use the landmark data to obtain more accurate estimates. 
This implies that we know approximately where, on the image, landmarks should appear- 
perhaps to within a few dozen picture elements (pixels). We shall ignore here any errors 
in our estimation of the orientations and sizes of the features; it is reasonable to assume 
that if the error in estimated position is small, then the errors in estimated orientation and 
size are negligible. 

If we regard the orientations and sizes of the landmarks on the image as shown, and their 
positions as approximately known, then the problem of landmark identification and loca- 
tion can be solved by a template matching process. This basically involves comparing a 
template or small image of the landmark with the picture, in the range of positions where 
the landmark could be located. If a good match is obtained, the landmark has been detected, 
and the position of this good match gives the location of the landmark. 


142 



There are many different template matching processes that could be used for this purpose, 
and the matching can be implemented in various ways (digitally, optically, or by a human 
operator). We shall consider here only automated digital techniques. 

The objective of this proposed effort is to make the template matching as computationally 
efficient as possible, while keeping machine storage requirements at the minicompuier level. 

The standard digital template matching technique, which is based on cross-correlation, is 
relatively slow. Substantially faster results can be obtained by using an absolute difference 
measure of mismatch (sometimes known as the sequential similarity detection algorithm 
(SSDA). This approach is not only cheaper computationally, but also lends itself to the 
u^ of efficient search techniques for speeding up the detection of mismatches, as shown 
by Nagel and Rosenfeld, 1972. It can be combined with various pre- and post-processing 
techniques that should increase the sharpness of the matches obtained, thus increasing the 
accuracy, as well as the speed, with which landmarks can be located. 

TEMPLATE MATCHING 


Use of the Absolute Difference Misnnatch Measure 

The most commonly used measure of the match between a template, t (x, y), and a oicture 
p (x, y), is the correlation coefficient (or normalized cross-correction) 



(x, y) p (x + u, y + v) dxdy 



(x, y) dxdy 



(x + u, y + v) dxdy 


where (u, v) is the displacement of the template relative to the picture. Here the area of 
integration, T, is the area occupied by the template. This coefficient takes on values be- 
tween 0 and 1 , and has value 1 only whe^^ t (x, y) and p (x + u, y + v) are identical (over 
the area T) except for a multiplicative constant. In the digital case, the expression for the 
correlation coefficient is identical, except that the integrals are replaced by sums over pic 
ture element arrays. 

Computation of the corr :!at*'on coefficient is relatively slow, because it involves a large 
number of multiplications of template values by picture values for each position of the 
template relative to the picture. The process can be speeded up by using Fourier transform 
techniques to compute the cross-correl?tion; however, this approach is more costly in com- 
puter memory requirements, and will not be discussed further here. 

A more quickly computable measure of the mismatch between a template and a picture is 
the integrated (or summed) absolute difference. 
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y>-p(x + u, y + v) 


dxdy 


This measure is zero wh 'n t -ind p are identical (over the area T), Computation of this 
measure involves no multiplications; many fewer arithmetic operations are required than 
in the case of the corre’ation measure (for a quantitative comparison see Bainea and Silver- 
man, 1972). In addition, use of the absolute difference measure makes it possible to take 
advantage of significant shortcuts in the matching process. 


Stoppiiitt Criteria 

A major advantage of the sum of the absolute dii fcrences as a mismatch measure is that 
this measure grows monotonically with the number of points whose absolute differences 
have been added into the sum. This means that if we are looking for a point of best match 
in j given region, we can stop adding up absolute differences in a given position (u, v) as 
soon as the sum has grown larger than the smallest sum sc far obtained in any position 
(u^, ). since we now know that (u, v) ernnot possibly be the position of best match. 

We can do even better if we set up an absolute tltreshold, such that if the mismatch measure 
exceeds & at a given point we will not accept that point as a pouit of match. It is now 
sufficient to add up terms of the sum until 6 is exceeded : when this happens, we can reject 
the given point and move to ^he next point. 

By using both mismatch comparison and absolute thresholding, the process of finding a 
best match (below .he mismatch threshold) can be speeded up considerably. For points 
where the match is poor, the threshold should be exceeded rapidly so that oniy a fraction 
of me terms in the sum need be computed for such points. Even for points where the match 
is good, we need not always compute all the terms of the sum, since we can slop as soon as 
we exceed the smallest previously obtained sum. 


Optimum Matching Sequences 

J( has been pointed out by Nagel and Rosenfeld (1972) that the expected amount of time 
required to exceed a mismatch threshold can be reduced by matching template pixels against 
the corresponding image pixels in a preselected order. To »Hustrate this idea, let us consider 
;i simple example. Suppose that the template and image contain only three gray levels, e.g., 
0, I , and 2, and that, in the given r^ion of tne image, these levels occur with relative fre^ 
quencies of 1 /2, 1 0, and 1 /6, respectively. If • mplate is not at or near i. * correct 
match position, we can assume that the :mage el that coincides with a given point 

of the template is uncorrelated with the template jZky level at that point. Thus, we can 
compute the expected absolute gray level difference between template and image at a point 
as follows: 
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Template Probability of 


Gray Level 


Difrerence 



0 


-1 

0 

1 



1/6 

1/3 

1/3 

0 

0 

1 

0 

1/6 

1/3 

1/2 

0 


0 

0 

1/6 

1/3 

•/2 


Probability of Absolute 
Difference 

Expected Absolute 
Difference 

0 1 



1/2 1/3 

1/6 

2^3 

1/3 2/3 

0 

213 

1/6 1/3 

1/2 

4/3 


In other words, the expected absolute difference between template and image at a template 
point having gray level 0 or 1 is 2/3; but. at a template point of gray level 2, the expected 
absolute difference is 4/3. Thus we can expect, on the average, a faster rate of growth of 
the sum of absolute differences if we compute the difference first for template points that 
have gray level 2. Thus, the mismatch threshold (or the lowest level of mismatch previously 
foun"*) is likely to be exceeded sooner if we compare template points with image points in 
a special order-namely. first using template points of gray level 2 -rather than using the 
template points in an arbitrary order* 

( meraiizing the example just given, it is easily seen that the mismatch threshold (absolute 
or relative) will be exceeded faster on the average if the template points are compared with 
the corresponding image points in decreasing order of expected absolute gray level differ- 
ence. The degree of speed-up achievable in this way depends, of course, on the probability 
distribution of gray levels in ;he given region of the image: if all gray levels are equally 
likely, no speedup is possible by this method. 

Tlie method just described should be even more effective when we are matching midti- 
spectral, rather than gray scale, images. This is because a probability distribution of multi- 
spectral vectc’* values should be even more nonuniform than a distribution of gray levels 
would be. Thus it should be possible to select template points whose expected absolute 
differences from the corresponding image points should be oMite large. Here the dift<^r- 
ences must, of course, be defined in vector icrms, e.g., as sums cf absolute differences of 
the individual vector components. 

The speedup obtainable by this method will not be the same for all landmarks. But the 
best results should be obtained for landmarks in whose vicinity the distribution of gray 
levels, or multispectral values, is higltly nonuniform. Given a set of available landmarks, 
one can determine, for each of them, how much speedup can be expected, and one can 
select preferred landmarks for which the expected speedup is greatest. 


* It may take more time to locate th». Jesired template points than it would if we tested them, perhaps in a simple 
raster sequence; but Nagel and Rosenfeld (1972) have shown that their method is faster than raster matching in 
spite of this. 
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The position of best match, determined by any method, may not be sharply defined, since 
there may be a set of nearby positions for which the match is nearly as good. It lias been 
known for at least 20 years that match sharpness can be improved by differentiating or 
by high-pass Filtering the image and the template before matching them. In other words, 
outlines of regions give sharper matches than do solid regions.* It would be of interest to 
investigate whether the matching speedup process described above is improved by using 
differentiated or outlined images and templates. This will certainly be the case if the 
differentiated image has a less uniform distribution of gray levels (or spectral values) than 
the original image, as we would normally expect. 

ITERATIVE RESAMPLING 

Match positions can be determined to within closer than a pixel by resampling the image 
(at a grid of points that have fractional coordinates, :n terms of the original sampling grid) 
and assigning gray levels to the new pixels by interpolation. The match at such an inter- 
polated position may be better than the match at any of the original (intercoordinate) 
position. The investigation of interpolated match positions is probably best done as a fine 
tuning of a coarse match position found without use of interpolation. The interpolated 
images will tend to have more uniform gray level or spectral distributions than the original 
image, so that the matching speedup process will probably not be as effective during this 
fine tuning stage as at the coarse search stage. 

A refinement allowing more accurate (part pixel) registration, when such accuracy is 
warranted, may be described as follows: 

Let F (x, y) be the cross-correlation function in terms of picture displacements (x, y) from 
the assumed best match point, x = 0, y = 0. Fit a quadratic to F (x, y) of the form 

f y) “ ax^ + 2h + by^ + 2gx + 2fy 

We may determine the coefficients from the values of F at a few points. At x = n, y = 0, 
for example; 


Similarly, 


F(n, 0) = an^ + 2gn 


F (-n, 0) - an^ - 2gn 


hence 


a = { F (n, 0) + F (-n, 0)]/2n^ 
g = ( F (n, 0) - F (-n, 0)] /4n 


* It may be of further advantage to use a thresholded outline image so as to reduce the effects of overall gray level 
differences due to seasonal changes, and so forth. 
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We can also find 


b = [ F (0, n) + F (0. -n)] 2n* 
r = I F(0, n)-F(0,-n)l/4n 
h = [ F(n, n) - F(n. 0) - F (0, n)l/2n^ 

Differentiating. F (x, y) is a maximum when* 

3^0 * hy„ + g = 0 

and 

so that, 

x„= (hf-bgV(ab-h^) 
y^= (gh -af)/(ab -h^) 

Then shift origin to (x^. y^) and repeat with n/2 area size. Since only 5 values of F need to 
be obtained per iteration, the speed of such a procedure may prove to be adequate. 

As an example in one dimension, consider the following template: 

3 19 13 

and the image line gray values: 

6 6 4 8 5 6 6 

Take (F (x)) = ax^ + 2gx. 

(3X 6) + (l X 6)t(9X 4)-t-(l X 8) + (3X 5) 

■v/(3^ + 1 ^ + 9^ + 1 ^ + 3^) (6^ + 6^ + 4^ + 8^ + 5^) 

= 0.6208 = a-2g 

(3 X 4) + (l X 8) + (9X 5) + (l X 6)+(3 X 6) 

F(l) = — - — — . - - 

>/( 3 ^ + 1 ^ + 9 ^ + 1 ^ + 3 ^) ( 4 ^ + 8 ^ + 5 ^ + 6 ^ + 6 ^) 

= 0.6656 = a + 2g 
a = 1/2(F(1) + F(-1)) = 0.6432 
g = 1/4(F(1)- F(-l)) = 0.0112 


After determining coefficients, we require that F be positive definite, i.e., ab > h^. to guarantee a true maximum of F. 
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Differentiating. F (x) is a maximum when ax^ + g = 0, i.e., for x^ = -g/a = -0.01 74 we see; 


Gray Scale 

Old Coordinate 

New Coordinate 

6 

-3 

-2.9826 

6 


-1.9826 

4 

-1 

-0.9826 

8 

0 

0.0174 

5 

1 

1.0174 

6 


2.0174 

6 

3 

3.0174 

In this new coordinate system we will need F(-0.5), hence 
the gray level values at X = -2.5. -1 .5. -0.5, 0.5. 1 .5, 2.5. 

by interpolation we will need 


If linear interpolation is used, then 


(G(Xj)-G(x,)) 

G(x)= G(x )+ ; — (x-x ) 

(Xj-X,) 

thus. 


G(-2.5)= 6.0000 
G(-I.5)= 5.0348 
G(-0.5)= 5.9304 
G ( 0.5) = 6.5522 


and so forth. 

Then a new iteration is undertaken and the process continues until the origin shift indicated 
is satisfactorily small. 

OTHER SPEEDUPS 

A possibility for still further improvement in matching efficiency, not discussed here up to 
now. is to identify locations in the image where the match can be expected to be good. This 
is typically done by making some simple local measurements on the image, and finding 
positions where the values of these measurements are close to their values for the template. 

In the present application, the positional uncertainty of the landmarks is not expected to 
be very large, so that it seems hkely that much will be gained by the use of this approach. 

There is. however, an important possibility for reducing the number of match positions that 
need be tried, once a match has been found f jr the first landmark. If this match is correct, 
the positional certainty of the remaining landmarks has now been reduced, since there 
are now fewer degrees of freedom. This is true even if the match is regarded as only approxi- 
mate. Thus, after the first match is found, the next match can be searched for over a 
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ivlaii\ol\ smaller region. Alter enough matches have been toiind. additional matches liave 
vcr\ little positional uncertainty, so tliat we can search for them over very small regions. 

Of course, if the subsequent matches are not found in the expected places, they must be 
sought for over wider regions, and tlie original matches must be reevaluated. Tints a feed- 
back process can be used to zero in on a best combination of match positions. 

TECHNICAL DISCUSSION 

The specific application to which tnis rcscarcli is directed is the navigation of the Synchronous 
Meteorological Satellite (SMS I. In umc ii is lioped that the processing of recognized land- 
marks for navigation (both for location error correction and relocation) can take place in 
quasi-real-time, which in this context is approximately 5 minutes. Before this can come 
about, there is a need to determine landmark processing algorithms that are rapid, accurate, 
reliable, and economical in storage. 

It is recognized that several operational questions influence tlie utility of the algorithm 
development results. In particular, eventual implementation may be effected on a mini- 
computer with a hardware dot product. The availability of such a feature would have 
some significance to the efficiency of algorithms in the multispectral domain. 

Also, error budgeting will play a central role. It may not be useful, for example, to process 
the image data with extreme accuracy since the errors in the landmark surveys may well 
predominate (current landmark accuracies are believed to be about 1/iOOth of a degree in 
latitude and longitude, and may have a systematic bias). Hence it may be useful to con- 
centrate on using the same landmarks repetitively for relative image -to-i mage geometric 
transformations rather than absolute geodetic location determination. 

Thirdly, cloudcover distribution probability functions will affect the ability to use the land- 
marks. The geometric relationships require a well distributed set of control points if orien- 
tation and orbit parameters are to be determined accurately. Inadequate or ill-distributed 
control point data give rise to ill-conditioned matrices and the calculations become unreliable. 
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CLOSED FORM SATELUTE THEORY WITH 
EQUINOCTIAL ELEMENTS 


R. Broucke 

University of Texas 
Austin. Texas 


A computer software system has been developed to generate the closed-form literal first- 
order perturbations due to any harmonic in the potential. In a first approach, classical 
elements and the true anomaly are used. In another approach, equinoctial elements and the 
true longitude are the basic variables. The solution in equinoctial elements does not have 
any zero eccentricity or zero ir 'ination singularities. 

In the case of tesseral and zonal harmonics, the rotation of the central body is neglected. 
The expansion of the potential is done with simple recurrence formulas. 
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FORMULATION OF AN ARBITRARY GEOPOTENTIAL TERM 
(TESSERAL) IN EQUINOCTIAL VARIABLES 


P, Cefola 

Computer Sciences Corporation 
Silver Spring, Maryland 


Previously, general formulas for the averaged disturbing potential were obtained ( AlAA pre- 
print 75-9) in equinoctial elements for the zonal and third-body harmonics, in addition, 
methods were given for recursive computation of these potentials. The current paper extends 
the applicability of this model by giving an explicit expression in equinoctial elements for 
the arbitrary geopotential term: 



P (sin 6 ) (C cos m X + S _sin m X] 

nm ‘ nm nm * 


Expressions for the spherical harmonics P^ ^ (sin 0) cos mX and P^ ^ (sin <t>) sin mX are 
obtained in terms of Jacobi polynomials with the argument (1 - p^ -q^ )/(l + p^ + q^ ) 

(p and q are equinoctial elements), polynomial functions which are straightforward gen- 
eralizations of the and polynomials appearing in the zonal potential, the true longi- 
tude, and the Greenwich sidereal time, 6 is fixed during the averaging period. However, 
an expansion of the true longitude in terms of the mean longitude and the equinoctial 
elements k and h would allow the consideration of resonant cases. Finally, consideration 
is given to recursive computation of the averaged potential for tesserals. 
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AN ANALYTIC METHOD TO ACCOUNT FOR ATMOSPHERIC 
DRAG WITH AN ANALYTIC SATELLITE THEORY 


N. Bonavito and R. Gordon 

Goddard Space Flight Center 
Greenbelt, Maryland 


The motion of an artificial earth satellite in the presence of air drag and the eartlTs gravita- 
tional potential is considered. In contrast to the classical methods of numerical integration, 
this approach presents a quadrature algorithm employing analytical expressions for the 
variation of orbital elements produced by the air drag. These expressions are w'ell-defined 
over expanded subintervals of the solution, and produce accurate agreement with profiles 
of tabular density. This procedure then allows a flexibiliiy in the selection of end points 
of the subintervals, which in turn ensures a minimum error bound on the required analytical 
function. 

In this method the effect of oblateness is accounted for by either the Vinti spheroidal theory, 
the Brouwer orbit theory, or the Brouwer-Lydanne theory. The changes due to atmospheric 
resistance for a nonrotating sphere are accounted for by the solutions of the variational 
equations, evaluated with the appropriate theory. 
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